### Replication Data and Code for: "Comment on: ### 
### Cropland abandonment in the context of drought, economic### 
###restructuring, and migration in Northeast Syria"### 
###by Konstantin Ash and Duy Trinh### 
###Forthcoming in Environmental Research Letters### 

#### Section 1: Setting up work environment ----
#### Set working Directory
setwd("~/ash_trinh_erl_comment_replication/")
~~

#### Install essential packages
install.packages("sf")
install.packages("gdal")
install.packages("tidyr")
install.packages("dplyr")
install.packages("raster")
install.packages("e1071")
install.packages("Rtools")

#### Load Libraries
library(sf)
library(rgdal)
library(tidyr)
library(dplyr)
library(raster)

#### Section 2: Loading Data ----

# Downloading and preparing the administrative boundaries file
# The Adm were obtained from the GADM database, managed and hosted by UC Davis.

# Adm.URL<- "https://geodata.ucdavis.edu/gadm/gadm4.1/shp/gadm41_SYR_shp.zip"
# dwnload<-"Input/adm.zip"
# Zip<-(download.file(Adm.URL,dwnload))
# unzip("Input/adm.zip",exdir = "Input")
#### Load administrative boundaries
adm.country<-st_read("~/ash_trinh_erl_comment_replication/gadm41_SYR_shp/gadm41_SYR_0.shp")

#### Load all raster images
# Fallowraster<-raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow2001.tif")
fallownew<-raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2000.tif")

Fallow2000<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2000.tif"))
Fallow2001<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2001.tif"))
Fallow2002<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2002.tif"))
Fallow2003<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2003.tif"))
Fallow2004<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2004.tif"))
Fallow2005<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2005.tif"))
Fallow2006<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2006.tif"))
Fallow2007<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2007.tif"))
Fallow2008<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2008.tif"))
Fallow2009<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2009.tif"))
Fallow2010<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2010.tif"))
Fallow2011<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2011.tif"))
Fallow2012<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2012.tif"))
Fallow2013<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2013.tif"))
Fallow2014<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2014.tif"))
Fallow2015<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2015.tif"))
Fallow2016<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2016.tif"))

Fallow_stack<-array(c(Fallow2000,Fallow2001,Fallow2002,Fallow2003,Fallow2004,Fallow2005,Fallow2006,Fallow2007,
                      Fallow2008,Fallow2009,Fallow2010,Fallow2011,Fallow2012,Fallow2013,Fallow2014,Fallow2015,
                      Fallow2016),
                    dim=c(5554,9762,17))

rm(Fallow2000,Fallow2001,Fallow2002,Fallow2003,Fallow2004,Fallow2005,Fallow2006,Fallow2007,Fallow2008,Fallow2009,Fallow2010,Fallow2011,Fallow2012,Fallow2013,Fallow2014,Fallow2015,Fallow2016)

Abandon<-matrix(NA,nrow = 5554,ncol=9762)
dim(Abandon)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster<- raster(Abandon)
crs(Abandon.raster)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster<-setExtent(Abandon.raster,fallownew)
writeRaster(Abandon.raster,"Abandoned_Final.tif", datatype='INT1U')

Abandon.raster %>%
  as.vector() %>%
  table() 

##Isolate Values###

Abandon.raster[Abandon.raster>3]<-NA

writeRaster(Abandon.raster,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_Final0002.tif", datatype='INT1U')

Abandon.raster<- raster(Abandon)

plot(fallownew)

###2000-2002

fallownew<-raster("NewFallow/Fallow_2000.tif")

Fallow2000<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2000.tif"))
Fallow2001<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2001.tif"))
Fallow2002<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2002.tif"))

Fallow_stack<-array(c(Fallow2000,Fallow2001,Fallow2002
),
dim=c(5554,9762,17))

rm(Fallow2000,Fallow2001,Fallow2002)


Abandon0002<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0002)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0002[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0002[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0002<- raster(Abandon0002)
crs(Abandon.raster0002)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0002<-setExtent(Abandon.raster0002,fallownew)
writeRaster(Abandon.raster0002,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0002.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0002 %>%
  as.vector() %>%
  table() 

plot(Abandon.raster0002)

plot(fallownew)
agriculture<-fallownew
agriculture[agriculture>0]<-0
plot(agriculture)

agriculture[Abandon.raster0002 == 3] <- 1
Abandon0002<-agriculture


###2001-2003

Fallow2001<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2001.tif"))
Fallow2002<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2002.tif"))
Fallow2003<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2003.tif"))

Fallow_stack<-array(c(Fallow2001,Fallow2002,Fallow2003
),
dim=c(5554,9762,17))

rm(Fallow2001,Fallow2002,Fallow2003)


Abandon0103<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0103)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0103[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0103[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0103<- raster(Abandon0103)
crs(Abandon.raster0103)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0103<-setExtent(Abandon.raster0103,fallownew)
writeRaster(Abandon.raster0103,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0103.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0103 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0103 == 3] <- 1
Abandon0103<-agriculture

#2002-2004

Fallow2002<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2002.tif"))
Fallow2003<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2003.tif"))
Fallow2004<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2004.tif"))

Fallow_stack<-array(c(Fallow2002,Fallow2003,Fallow2004
),
dim=c(5554,9762,17))

rm(Fallow2002,Fallow2003,Fallow2004)


Abandon0204<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0204)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0204[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0204[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0204<- raster(Abandon0204)
crs(Abandon.raster0204)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0204<-setExtent(Abandon.raster0204,fallownew)
writeRaster(Abandon.raster0204,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0204.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0204 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0204 == 3] <- 1
Abandon0204<-agriculture



#2003-2005

Fallow2003<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2003.tif"))
Fallow2004<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2004.tif"))
Fallow2005<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2005.tif"))

Fallow_stack<-array(c(Fallow2003,Fallow2004,Fallow2005
),
dim=c(5554,9762,17))

rm(Fallow2003,Fallow2004,Fallow2005)


Abandon0305<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0305)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0305[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0305[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0305<- raster(Abandon0305)
crs(Abandon.raster0305)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0305<-setExtent(Abandon.raster0305,fallownew)
writeRaster(Abandon.raster0305,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0305.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0305 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0305 == 3] <- 1
Abandon0305<-agriculture


#2004-06

rm(Abandon.raster0002,Abandon.raster0103,Abandon.raster0204,Abandon.raster0305)

Fallow2004<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2004.tif"))
Fallow2005<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2005.tif"))
Fallow2006<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2006.tif"))

Fallow_stack<-array(c(Fallow2004,Fallow2005,Fallow2006
),
dim=c(5554,9762,17))

rm(Fallow2004,Fallow2005,Fallow2006)


Abandon0406<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0406)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0406[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0406[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0406<- raster(Abandon0406)
crs(Abandon.raster0406)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0406<-setExtent(Abandon.raster0406,fallownew)
writeRaster(Abandon.raster0406,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0406.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0406 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0406 == 3] <- 1
Abandon0406<-agriculture

rm(Abandon.raster0406)

#2005-2007

Fallow2005<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2005.tif"))
Fallow2006<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2006.tif"))
Fallow2007<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2007.tif"))

Fallow_stack<-array(c(Fallow2005,Fallow2006,Fallow2007
),
dim=c(5554,9762,17))

rm(Fallow2005,Fallow2006,Fallow2007)


Abandon0507<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0507)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0507[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0507[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0507<- raster(Abandon0507)
crs(Abandon.raster0507)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0507<-setExtent(Abandon.raster0507,fallownew)
writeRaster(Abandon.raster0507,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0507.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0507 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0507 == 3] <- 1
Abandon0507<-agriculture

rm(Abandon.raster0507)

#2005-2007

Fallow2005<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2005.tif"))
Fallow2006<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2006.tif"))
Fallow2007<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2007.tif"))

Fallow_stack<-array(c(Fallow2005,Fallow2006,Fallow2007
),
dim=c(5554,9762,17))

rm(Fallow2005,Fallow2006,Fallow2007)


Abandon0507<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0507)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0507[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0507[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0507<- raster(Abandon0507)
crs(Abandon.raster0507)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0507<-setExtent(Abandon.raster0507,fallownew)
writeRaster(Abandon.raster0507,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0507.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0507 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0507 == 3] <- 1
Abandon0507<-agriculture

rm(Abandon.raster0507)


#2006-2008

Fallow2006<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2006.tif"))
Fallow2007<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2007.tif"))
Fallow2008<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2008.tif"))

Fallow_stack<-array(c(Fallow2006,Fallow2007,Fallow2008
),
dim=c(5554,9762,17))

rm(Fallow2006,Fallow2007,Fallow2008)


Abandon0608<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0608)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0608[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0608[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0608<- raster(Abandon0608)
crs(Abandon.raster0608)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0608<-setExtent(Abandon.raster0608,fallownew)
writeRaster(Abandon.raster0608,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0608.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0608 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0608 == 3] <- 1
Abandon0608<-agriculture

rm(Abandon.raster0608)


#2007-2009

Fallow2007<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2007.tif"))
Fallow2008<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2008.tif"))
Fallow2009<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2009.tif"))

Fallow_stack<-array(c(Fallow2007,Fallow2008,Fallow2009
),
dim=c(5554,9762,17))

rm(Fallow2007,Fallow2008,Fallow2009)


Abandon0709<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0709)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0709[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0709[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0709<- raster(Abandon0709)
crs(Abandon.raster0709)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0709<-setExtent(Abandon.raster0709,fallownew)
writeRaster(Abandon.raster0709,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0709.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0709 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0709 == 3] <- 1
Abandon0709<-agriculture

rm(Abandon.raster0709)


#2008-2010

Fallow2008<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2008.tif"))
Fallow2009<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2009.tif"))
Fallow2010<-as.matrix(raster("~/ash_trinh_erl_comment_replication/NewFallow/Fallow_2010.tif"))

Fallow_stack<-array(c(Fallow2008,Fallow2009,Fallow2010
),
dim=c(5554,9762,17))

rm(Fallow2008,Fallow2009,Fallow2010)


Abandon0810<-matrix(NA, nrow = 5554, ncol=9762)

dim(Abandon0810)
for(i in 1:5554){
  for(j in 1:9762){
    ts<- Fallow_stack[i,j,]
    fcount<-0
    for(t in 1:dim(Fallow_stack)[3]){
      if (is.na(Fallow_stack[i,j,t])){
        Abandon0810[i,j]<-NA
        break
      }else{
        if (Fallow_stack[i,j,t]==1){
          fcount<-fcount+1
        }else{
          fcount<-0
        }
        if (fcount==3){
          Abandon0810[i,j]<-t
          break
        }
      }
    }
  } 
}

Abandon.raster0810<- raster(Abandon0810)
crs(Abandon.raster0810)<- "+proj=longlat +datum=WGS84 +no_defs"
Abandon.raster0810<-setExtent(Abandon.raster0810,fallownew)
writeRaster(Abandon.raster0810,"~/ash_trinh_erl_comment_replication/Abandonment/Abandoned_0810.tif", datatype='INT1U',overwrite=TRUE)

Abandon.raster0810 %>%
  as.vector() %>%
  table() 

agriculture<-fallownew
agriculture[agriculture>0]<-0

agriculture[Abandon.raster0810 == 3] <- 1
Abandon0810<-agriculture

rm(Abandon.raster0810)


####Zonal Statistics####


library(sf)
##setwd("E:/OneDrive/Documents/GIS/Eklund et al 2022 Extension/Ash 2023 REC Replication/SYR_adm")
syr_adm3<-readOGR(dsn=".", layer="SYR_adm3_noflares")

abandon2002 <- extract(Abandon0002, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE) 
table2002 <- as.data.frame(abandon2002)
abandon2003 <- extract(Abandon0103, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE)
table2003 <- as.data.frame(abandon2003)
abandon2004 <- extract(Abandon0204, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE)
table2004 <- as.data.frame(abandon2004)
abandon2005 <- extract(Abandon0305, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE)
table2005 <- as.data.frame(abandon2005)
abandon2006 <- extract(Abandon0406, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE) 
table2006 <- as.data.frame(abandon2006)
abandon2007 <- extract(Abandon0507, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE)
table2007 <- as.data.frame(abandon2007)
abandon2008 <- extract(Abandon0608, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE)
table2008 <- as.data.frame(abandon2008)
abandon2009 <- extract(Abandon0709, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE) 
table2009 <- as.data.frame(abandon2009)
abandon2010 <- extract(Abandon0810, syr_adm3, fun=mean, small=TRUE, sp=TRUE, na.rm=TRUE) 
table2010 <- as.data.frame(abandon2010)

table2002<-table2002 %>% rename(
  abandon0002=Fallow_2000.Classification
)
table2002<-table2002[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0002")]

table2003<-table2003 %>% rename(
  abandon0103=Fallow_2000.Classification
)
table2003<-table2003[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0103")]

table2004<-table2004 %>% rename(
  abandon0204=layer
)
table2004<-table2004[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0204")]

table2005<-table2005 %>% rename(
  abandon0305=Fallow_2000.Classification
)
table2005<-table2005[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0305")]

table2006<-table2006 %>% rename(
  abandon0406=Fallow_2000.Classification
)
table2006<-table2006[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0406")]

table2007<-table2007 %>% rename(
  abandon0507=Fallow_2000.Classification
)
table2007<-table2007[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0507")]

table2008<-table2008 %>% rename(
  abandon0608=Fallow_2000.Classification
)
table2008<-table2008[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0608")]

table2009<-table2009 %>% rename(
  abandon0709=Fallow_2000.Classification
)
table2009<-table2009[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0709")]

table2010<-table2010 %>% rename(
  abandon0810=Fallow_2000.Classification
)
table2010<-table2010[c("NAME_EN","ADMIN2_EN","ADMIN1_EN","PCODE","abandon0810")]

abandon<-merge(table2002,table2003,by=c("PCODE","NAME_EN","ADMIN2_EN","ADMIN1_EN"))
abandon<-merge(abandon,table2004,by=c("PCODE","NAME_EN","ADMIN2_EN","ADMIN1_EN"))
abandon<-merge(abandon,table2005,by=c("PCODE","NAME_EN","ADMIN2_EN","ADMIN1_EN"))
abandon<-merge(abandon,table2006,by=c("PCODE","NAME_EN","ADMIN2_EN","ADMIN1_EN"))
abandon<-merge(abandon,table2007,by=c("PCODE","NAME_EN","ADMIN2_EN","ADMIN1_EN"))
abandon<-merge(abandon,table2008,by=c("PCODE","NAME_EN","ADMIN2_EN","ADMIN1_EN"))
abandon<-merge(abandon,table2009,by=c("PCODE","NAME_EN","ADMIN2_EN","ADMIN1_EN"))
abandon<-merge(abandon,table2010,by=c("PCODE","NAME_EN","ADMIN2_EN","ADMIN1_EN"))



write.csv(abandon, "~/ash_trinh_erl_comment_replication/abandon.csv")

####Section 2: With Completed Data, Analysis Can Start Here#####


library(foreign)
library(lfe)
library(data.table)
library(ggplot2)
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(modelr)
abandon_main<-read.dta("~/ash_trinh_erl_comment_replication/main_abandon.dta")


abandon_main <- abandon_main[!grepl("Arbin|Damascus|Hajar Aswad|Jaramana", abandon_main$name_en), ]

setDT(abandon_main)
setkey(abandon_main, name_en, year)
abandon_main[, prcp.dev := shift(prcp - meanprcp), by=.(name_en)]
abandon_main[, temp.dev := shift(temp - meantemp), by=.(name_en)]



##Fix Stargazer##

detach("package:stargazer",unload=T)
# Delete it
remove.packages("stargazer")
# Download the source
download.file("https://cran.r-project.org/src/contrib/stargazer_5.2.3.tar.gz", destfile = "stargazer_5.2.3.tar.gz")
# Unpack
untar("stargazer_5.2.3.tar.gz")
# Read the sourcefile with .inside.bracket fun
stargazer_src <- readLines("stargazer/R/stargazer-internal.R")
# Move the length check 5 lines up so it precedes is.na(.)
stargazer_src[1990] <- stargazer_src[1995]
stargazer_src[1995] <- ""
# Save back
writeLines(stargazer_src, con="stargazer/R/stargazer-internal.R")
# Compile and install the patched package
install.packages("stargazer", repos = NULL, type="source")

library(stargazer)


continuousinteraction <- function(felmmodel, lmmodel, rhs1, intxn, minimum="min", maximum="max", miny=-8, maxy=1, incr="default", num_points = 1000, conf=.95, title="Interaction marginal effects plot", xlabel="Value of Intxn", ylabel="Estimated Marginal Effect", maincolor='black', secondcolor='black'){
  
  # Get coefficients
  beta_1 = felmmodel$coefficients[row.names(felmmodel$coefficients)==rhs1]
  
  # Get coefficient if interaction term==100
  beta_3 = felmmodel$coefficients[row.names(felmmodel$coefficients)==paste0(rhs1, ':', intxn)]
  
  # Extract Variance Covariance matrix
  covMat = felmmodel$clustervcv
  
  # Extract the data frame of the model
  mod_frame = model.frame(lmmodel)
  
  # Set range of the interaction variable
  # Minimum
  if (minimum == "min"){
    min_val = min(mod_frame[[intxn]])
  }else{
    min_val = minimum
  }
  # Maximum
  if (maximum == "max"){
    max_val = max(mod_frame[[intxn]])
  }else{
    max_val = maximum
  }
  
  # Check if minimum smaller than maximum
  if (min_val > max_val){
    stop("Error: Minimum X value greater than maximum value.")
  }
  
  # Determine intervals between values of the x variable
  if (incr == "default"){
    increment = (max_val - min_val)/(num_points - 1)
  }else{
    increment = incr
  }
  
  # Create list of intxn values at which marginal effect is evaluated
  intxn_2 <- seq(from=min_val, to=max_val, by=increment)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_3*intxn_2
  
  # Compute variances
  var_1 = covMat[rhs1,rhs1] + (intxn_2^2)*covMat[paste0(rhs1, ':', intxn), paste0(rhs1, ':', intxn)] + 2*intxn_2*covMat[rhs1,  paste0(rhs1, ':', intxn)]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  d <- as.data.table(cbind(intxn_2, delta_1, upper_bound, lower_bound))
  int <- data.frame(inter = mod_frame[[intxn]])
  
  #Plot
  ggplot(environment = environment()) +
    ggtitle(title) +
    geom_line(data=d, aes(y=delta_1, x=intxn_2), lwd=1.5, color='black') +
    geom_line(data=d, aes(y=upper_bound, x=intxn_2), lty=2, lwd=1.5, color=maincolor) +
    geom_line(data=d, aes(y=lower_bound, x=intxn_2), lty=2, lwd=1.5, color=maincolor) +
    geom_hline(yintercept=0, lty=3) +
    #geom_rug(data=int, mapping = aes(x=inter), alpha=.05, color=secondcolor) +
    ylab(ylabel) +
    xlab(xlabel) +
    coord_cartesian(ylim = c(miny,maxy), xlim=c(min_val,max_val)) +
    theme(title = element_text(size=16, angle=0, vjust=.5, face="bold"),
          plot.title = element_text(hjust=0.5),
          axis.title.x = element_text(size=18, angle=0, vjust=-0, face="bold"),
          axis.title.y = element_text(size=18,angle=90, vjust=1, face="bold"),
          axis.text.y = element_text(size=16,angle=0, vjust=.5, face="bold") ,
          axis.text.x = element_text(size=16, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(),
          legend.title=element_blank()
    )
  
}


###Table 1: Models for Effect of Abandonment on Nighttime Lights
drought.abandon.felm <- felm(lights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2008  & abandon_main$year <= 2010,], exactDOF=TRUE)
drought.abandon.somefe <- felm(lights ~ abandon_pct | admin1_en + year | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2008 & abandon_main$year <= 2010,])
predrought.abandon.felm <- felm(lights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year < 2008,], exactDOF=TRUE)
predrought.abandon.somefe <- felm(lights ~ abandon_pct | admin1_en + year | 0 | admin2_en, data=abandon_main[abandon_main$year < 2008,])
all.abandon.felm <- felm(lights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main, exactDOF=TRUE)
all.abandon.somefe <- felm(lights ~ abandon_pct | admin1_en + year | 0 | admin2_en, data=abandon_main)

summary(drought.abandon.felm)
summary(drought.abandon.somefe)
summary(predrought.abandon.felm)
summary(predrought.abandon.somefe)
summary(all.abandon.felm)
summary(all.abandon.somefe)


stargazer(all.abandon.felm, predrought.abandon.felm, drought.abandon.felm,
          type='latex',
          title = 'Land Abandonment and Nighttime Lights',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Nighttime Lights',
          covariate.labels = c( 'Land Abandonment Pct.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)

##### Table 2: Effect of Lagged Nighttime Lights on Abandonment####

abandon_main[, lights_lag := shift(lights), by=.(name_en)]
abandon_main[, lights_lag2 := shift(lights_lag), by=.(name_en)]


drought.abandon.felm.rev <- felm(abandon_pct ~ lights+lights_lag | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2008  & abandon_main$year <= 2010,], exactDOF=TRUE)
predrought.abandon.felm.rev <- felm(abandon_pct ~ lights+lights_lag | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year < 2008,], exactDOF=TRUE)
all.abandon.felm.rev  <- felm(abandon_pct~ lights+lights_lag | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main, exactDOF=TRUE)

summary(drought.abandon.felm.rev)
summary(predrought.abandon.felm.rev)
summary(all.abandon.felm.rev)

stargazer(all.abandon.felm.rev, predrought.abandon.felm.rev, drought.abandon.felm.rev,
          type='latex',
          title = 'Lagged Nighttime Lights and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Land Abandonment Pct.',
          covariate.labels = c( 'Nighttime Lights', 'Nighttime Lights in prior year'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)






###############################################################################
###Appendix 1: Operationalizing Land Abandonment based on First Year Fallowness    
###############################################################################

##### Using first year of Abandonment instead of third year####

abandon_main[, abandon_pct_secondyear := shift(abandon_pct, type="lead"), by=.(name_en)]
abandon_main[, abandon_pct_firstyear := shift(abandon_pct_secondyear, type="lead"), by=.(name_en)]



###Table A1: Climatic Stress and First Year Land Abandonment###

drought.felm_abandon_firstyear <- felm(abandon_pct_firstyear ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main[abandon_main$year >= 2006 & abandon_main$year <= 2008,], exactDOF=TRUE, psdef=FALSE)
drought.lm_abandon_firstyear <- lm(abandon_pct_firstyear ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=abandon_main[abandon_main$year >= 2006 & abandon_main$year <= 2008,])
predrought.felm_abandon_firstyear <- felm(abandon_pct_firstyear ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main[abandon_main$year < 2006,], exactDOF=TRUE, psdef=FALSE)
predrought.lm_abandon_firstyear <- lm(abandon_pct_firstyear ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=abandon_main[abandon_main$year < 2006,])
all.felm_abandon_firstyear<- felm(abandon_pct_firstyear ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main, exactDOF=TRUE, psdef=FALSE)
all.lm_abandon_firstyear <- lm(abandon_pct_firstyear ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=abandon_main)

summary(drought.felm_abandon_firstyear)
summary(drought.lm_abandon_firstyear)
summary(predrought.felm_abandon_firstyear)
summary(predrought.lm_abandon_firstyear)
summary(all.felm_abandon_firstyear)
summary(all.lm_abandon_firstyear)

nobs(drought.felm_abandon_firstyear)
nobs(drought.lm_abandon_firstyear)
nobs(predrought.felm_abandon_firstyear)
nobs(predrought.lm_abandon_firstyear)
nobs(all.felm_abandon_firstyear)
nobs(all.lm_abandon_firstyear)

stargazer(all.felm_abandon_firstyear, predrought.felm_abandon_firstyear, drought.felm_abandon_firstyear,
          type='latex',
          title = 'Precipitation, Temperature, and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: First Year Abandonment Pct.',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','first-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)



drought <- continuousinteraction(felmmodel=drought.felm_abandon_firstyear,
                                 lmmodel=drought.lm_abandon_firstyear,
                                 rhs1="prcp.dev",
                                 intxn="temp.dev",
                                 minimum="min",
                                 maximum="max",
                                 miny=-0.025,
                                 maxy=0.025,
                                 incr="default",
                                 num_points = 1000,
                                 conf=.95,
                                 title="Drought (06/08-08/10)",
                                 xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                 ylabel="Effect of \u0394 Precipitation on Land Abandonment",
                                 maincolor='darkorange',
                                 secondcolor='darkgray')

predrought <- continuousinteraction(felmmodel=predrought.felm_abandon_firstyear,
                                    lmmodel=predrought.lm_abandon_firstyear,
                                    rhs1="prcp.dev",
                                    intxn="temp.dev",
                                    minimum="min",
                                    maximum="max",
                                    miny=-0.025,
                                    maxy=0.025,
                                    incr="default",
                                    num_points = 1000,
                                    conf=.95,
                                    title="Pre-Drought (00/02-05/07)",
                                    xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                    ylabel="Effect of \u0394 Precipitation on Land Abandonment",
                                    maincolor='darkorange',
                                    secondcolor='darkgray')

all <- continuousinteraction(felmmodel=all.felm_abandon_firstyear,
                             lmmodel=all.lm_abandon_firstyear,
                             rhs1="prcp.dev",
                             intxn="temp.dev",
                             minimum="min",
                             maximum="max",
                             miny=-0.025,
                             maxy=0.025,
                             incr="default",
                             num_points = 1000,
                             conf=.95,
                             title="Full Sample (00/02-08/10)",
                             xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                             ylabel="Effect of \u0394 Precipitation on Land Abandonment",
                             maincolor='darkorange',
                             secondcolor='darkgray')
## plot figure 1
cairo_pdf(file = 'landabandonment3.pdf', width = 13, height = 6.5)
grid.arrange(all, predrought, drought, ncol=3)
dev.off()



### Table A2: Models for Effect of Abandonment on Nighttime Lights

drought.abandon.felm_firstyear <- felm(lights ~ abandon_pct_firstyear  | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2006  & abandon_main$year <= 2008,], exactDOF=TRUE)
drought.abandon.somefe_firstyear  <- felm(lights ~ abandon_pct_firstyear  | admin1_en + year | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2006 & abandon_main$year <= 2008,])
predrought.abandon.felm_firstyear  <- felm(lights ~ abandon_pct_firstyear  | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year < 2006,], exactDOF=TRUE)
predrought.abandon.somefe_firstyear  <- felm(lights ~ abandon_pct_firstyear  | admin1_en + year | 0 | admin2_en, data=abandon_main[abandon_main$year < 2006,])
all.abandon.felm_firstyear  <- felm(lights ~ abandon_pct_firstyear  | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main, exactDOF=TRUE)
all.abandon.somefe_firstyear  <- felm(lights ~ abandon_pct_firstyear  | admin1_en + year | 0 | admin2_en, data=abandon_main)

summary(drought.abandon.felm_firstyear )
summary(drought.abandon.somefe_firstyear )
summary(predrought.abandon.felm_firstyear )
summary(predrought.abandon.somefe_firstyear )
summary(all.abandon.felm_firstyear )
summary(all.abandon.somefe_firstyear)

nobs(drought.abandon.felm_firstyear )
nobs(drought.abandon.somefe_firstyear )
nobs(predrought.abandon.felm_firstyear )
nobs(predrought.abandon.somefe_firstyear )
nobs(all.abandon.felm_firstyear )
nobs(all.abandon.somefe_firstyear)


stargazer(all.abandon.felm_firstyear, predrought.abandon.felm_firstyear, drought.abandon.felm_firstyear,
          type='latex',
          title = 'Land Abandonment and Nighttime Lights',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Last Year Nighttime Lights',
          covariate.labels = c( 'First Year Abandonment Pct.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)

##### Table A3: Effect of Lagged Nighttime Lights on Abandonment####




drought.abandon.felm.rev_firstyear <- felm(abandon_pct_firstyear ~ lights+lights_lag | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2006  & abandon_main$year <= 2008,], exactDOF=TRUE)
predrought.abandon.felm.rev_firstyear <- felm(abandon_pct_firstyear ~ lights+lights_lag | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year < 2006,], exactDOF=TRUE)
all.abandon.felm.rev_firstyear  <- felm(abandon_pct_firstyear~ lights+lights_lag | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main, exactDOF=TRUE)

summary(drought.abandon.felm.rev_firstyear)
summary(predrought.abandon.felm.rev_firstyear)
summary(all.abandon.felm.rev_firstyear)

nobs(drought.abandon.felm.rev_firstyear)
nobs(predrought.abandon.felm.rev_firstyear)
nobs(all.abandon.felm.rev_firstyear)

stargazer(all.abandon.felm.rev_firstyear, predrought.abandon.felm.rev_firstyear, drought.abandon.felm.rev_firstyear,
          type='latex',
          title = 'Lagged Nighttime Lights and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: First Year Abandonment Pct.',
          covariate.labels = c( 'Last Year Nighttime Lights', 'Second Year Nighttime Lights'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)

###############################################################################
###Appendix 2: Climatic Stress and Land Abandonment ###########################     
###############################################################################


###Table A4: Time by location fixed effects of temperature and precipitation on abandonment#####

drought.felm_abandon <- felm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main[abandon_main$year >= 2008 & abandon_main$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
drought.lm_abandon <- lm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=abandon_main[abandon_main$year >= 2008 & abandon_main$year <= 2010,])
predrought.felm_abandon <- felm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main[abandon_main$year < 2008,], exactDOF=TRUE, psdef=FALSE)
predrought.lm_abandon <- lm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=abandon_main[abandon_main$year < 2008,])
all.felm_abandon <- felm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main, exactDOF=TRUE, psdef=FALSE)
all.lm_abandon <- lm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=abandon_main)

summary(drought.felm_abandon)
summary(drought.lm_abandon)
summary(predrought.felm_abandon)
summary(predrought.lm_abandon)
summary(all.felm_abandon)
summary(all.lm_abandon)


stargazer(all.felm_abandon, predrought.felm_abandon, drought.felm_abandon,
          type='latex',
          title = 'Precipitation, Temperature, and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Land Abandonment Pct.',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','first-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)



####Table A5/Figure A1: Effect of Climatic Stress and Land Abandonment on Nighttime Lights####

drought.felm_full <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+abandon_pct | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main[abandon_main$year >= 2008 & abandon_main$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
drought.lm_full <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+abandon_pct, data=abandon_main[abandon_main$year >= 2008 & abandon_main$year <= 2010,])
predrought.felm_full <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+abandon_pct | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main[abandon_main$year < 2008,], exactDOF=TRUE, psdef=FALSE)
predrought.lm_full <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+abandon_pct, data=abandon_main[abandon_main$year < 2008,])
all.felm_full <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=abandon_main, exactDOF=TRUE, psdef=FALSE)
all.lm_full <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+abandon_pct, data=abandon_main)

summary(drought.felm_full)
summary(drought.lm_full)
summary(predrought.felm_full)
summary(predrought.lm_full)
summary(all.felm_full)
summary(all.lm_full)

stargazer(all.felm_full, predrought.felm_full, drought.felm_full,
          type='latex',
          title = 'Climatic Stress and Nighttime Lights',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Nighttime Lights',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', 'Abandonment Pct.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are multiway clustered on','second-level administrative districts and years.',
                    'Only the 185 nawahi Eklund et. al. 2024', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)


drought <- continuousinteraction(felmmodel=drought.felm_full,
                                 lmmodel=drought.lm_full,
                                 rhs1="prcp.dev",
                                 intxn="temp.dev",
                                 minimum="min",
                                 maximum="max",
                                 miny=-0.5,
                                 maxy=0.5,
                                 incr="default",
                                 num_points = 1000,
                                 conf=.95,
                                 title="Drought (06/08-08/10)",
                                 xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                 ylabel="Effect of \u0394 Precipitation on Nighttime Lights",
                                 maincolor='darkorange',
                                 secondcolor='darkgray')

predrought <- continuousinteraction(felmmodel=predrought.felm_full,
                                    lmmodel=predrought.lm_full,
                                    rhs1="prcp.dev",
                                    intxn="temp.dev",
                                    minimum="min",
                                    maximum="max",
                                    miny=-0.5,
                                    maxy=0.5,
                                    incr="default",
                                    num_points = 1000,
                                    conf=.95,
                                    title="Pre-Drought (00/02-05/07)",
                                    xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                    ylabel="Effect of \u0394 Precipitation on Nighttime Lights",
                                    maincolor='darkorange',
                                    secondcolor='darkgray')

all <- continuousinteraction(felmmodel=all.felm_full,
                             lmmodel=all.lm_full,
                             rhs1="prcp.dev",
                             intxn="temp.dev",
                             minimum="min",
                             maximum="max",
                             miny=-0.5,
                             maxy=0.5,
                             incr="default",
                             num_points = 1000,
                             conf=.95,
                             title="Full Sample (00/02-08/10)",
                             xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                             ylabel="Effect of \u0394 Precipitation on Nighttime Lights",
                             maincolor='darkorange',
                             secondcolor='darkgray')
## plot figure A1
cairo_pdf(file = 'landabandonment1.pdf', width = 13, height = 6.5)
grid.arrange(all, predrought, drought, ncol=3)
dev.off()



#### Table A6/Figure 2: Third Level Admin and Year Fixed Effects of Lagged Drought on Abandonment####



drought.somefelm_abandon <- felm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en + year  | 0 | admin1_en, data=abandon_main[abandon_main$year >= 2008 & abandon_main$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
predrought.somefelm_abandon <- felm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en + year | 0 | admin1_en, data=abandon_main[abandon_main$year < 2008,], exactDOF=TRUE, psdef=FALSE)
all.somefelm_abandon <- felm(abandon_pct ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en + year | 0 | admin1_en, data=abandon_main, exactDOF=TRUE, psdef=FALSE)

summary(drought.somefelm_abandon)
summary(predrought.somefelm_abandon)
summary(all.somefelm_abandon)



stargazer(all.somefelm_abandon, predrought.somefelm_abandon, drought.somefelm_abandon,
          type='latex',
          title = 'Precipitation, Temperature, and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Land Abandonment Pct.',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "No", "No", "No")),
          notes = c('Standard errors are clustered on','first-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)


drought <- continuousinteraction(felmmodel=drought.somefelm_abandon,
                                 lmmodel=drought.lm_abandon,
                                 rhs1="prcp.dev",
                                 intxn="temp.dev",
                                 minimum="min",
                                 maximum="max",
                                 miny=-0.025,
                                 maxy=0.025,
                                 incr="default",
                                 num_points = 1000,
                                 conf=.95,
                                 title="Drought (06/08-08/10)",
                                 xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                 ylabel="Effect of \u0394 Precipitation on Land Abandonment",
                                 maincolor='darkorange',
                                 secondcolor='darkgray')

predrought <- continuousinteraction(felmmodel=predrought.somefelm_abandon,
                                    lmmodel=predrought.lm_abandon,
                                    rhs1="prcp.dev",
                                    intxn="temp.dev",
                                    minimum="min",
                                    maximum="max",
                                    miny=-0.025,
                                    maxy=0.025,
                                    incr="default",
                                    num_points = 1000,
                                    conf=.95,
                                    title="Pre-Drought (00/02-05/07)",
                                    xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                    ylabel="Effect of \u0394 Precipitation on Land Abandonment",
                                    maincolor='darkorange',
                                    secondcolor='darkgray')

all <- continuousinteraction(felmmodel=all.somefelm_abandon,
                             lmmodel=all.lm_abandon,
                             rhs1="prcp.dev",
                             intxn="temp.dev",
                             minimum="min",
                             maximum="max",
                             miny=-0.025,
                             maxy=0.025,
                             incr="default",
                             num_points = 1000,
                             conf=.95,
                             title="Full Sample (00/02-08/10)",
                             xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                             ylabel="Effect of \u0394 Precipitation on Land Abandonment",
                             maincolor='darkorange',
                             secondcolor='darkgray')
## plot figure 2
cairo_pdf(file = 'landabandonment2.pdf', width = 13, height = 6.5)
grid.arrange(all, predrought, drought, ncol=3)
dev.off()




###############################################################################
###Appendix 3: MODIS Land Cover Percentage as Predictor of Land Abandonment####
###############################################################################

#setwd("~/ash_trinh_erl_comment_replication/MODIS")
library(terra)
modis01 <- rast('modis01.tif')
modis02 <- rast('modis02.tif')
modis03 <- rast('modis03.tif')
modis04 <- rast('modis04.tif')
modis05 <- rast('modis05.tif')
modis06 <- rast('modis06.tif')
modis07 <- rast('modis07.tif')
modis08 <- rast('modis08.tif')
modis09 <- rast('modis09.tif')
modis10 <- rast('modis10.tif')

m <- c(0, 12.1, 0,  12.9, 13.1, 1,  13.9, 17.1, 0)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
modis01 <- classify(modis01, rclmat)
modis02 <- classify(modis02, rclmat)
modis03 <- classify(modis03, rclmat)
modis04 <- classify(modis04, rclmat)
modis05 <- classify(modis05, rclmat)
modis06 <- classify(modis06, rclmat)
modis07 <- classify(modis07, rclmat)
modis08 <- classify(modis08, rclmat)
modis09 <- classify(modis09, rclmat)
modis10 <- classify(modis10, rclmat)

writeRaster(modis01,"modis01_urban.tif", datatype='INT1U')
writeRaster(modis02,"modis02_urban.tif", datatype='INT1U')
writeRaster(modis03,"modis03_urban.tif", datatype='INT1U')
writeRaster(modis04,"modis04_urban.tif", datatype='INT1U')
writeRaster(modis05,"modis05_urban.tif", datatype='INT1U')
writeRaster(modis06,"modis06_urban.tif", datatype='INT1U')
writeRaster(modis07,"modis07_urban.tif", datatype='INT1U')
writeRaster(modis08,"modis08_urban.tif", datatype='INT1U')
writeRaster(modis09,"modis09_urban.tif", datatype='INT1U')
writeRaster(modis10,"modis10_urban.tif", datatype='INT1U')


modis0110 <- modis10-modis01
plot(modis0110)


#setwd("~/ash_trinh_erl_comment_replication/SYR_adm")
syr_adm3<-readOGR(dsn=".", layer="SYR_ADM3_clip")

test <- vect("SYR_ADM3_clip.shp")
plot(modis01)
plot(modis10)
plot(syr_adm3, add=T)

syr_adm3_name<-as.data.frame(syr_adm3@data$Name)

urban2001 <- terra::zonal(modis01, test, fun=mean, na.rm=TRUE) 
landcover01 <- as.data.frame(urban2001)
main_urban<-bind_cols(syr_adm3_name, landcover01)
urban2002 <- terra::zonal(modis02, test, fun=mean, na.rm=TRUE)  
landcover02 <- as.data.frame(urban2002)
main_urban<-bind_cols(main_urban, landcover02)
urban2003 <- terra::zonal(modis03, test, fun=mean, na.rm=TRUE)  
landcover03 <- as.data.frame(urban2003)
main_urban<-bind_cols(main_urban, landcover03)
urban2004 <- terra::zonal(modis04, test, fun=mean, na.rm=TRUE) 
landcover04 <- as.data.frame(urban2004)
main_urban<-bind_cols(main_urban, landcover04)
urban2005 <- terra::zonal(modis05, test, fun=mean, na.rm=TRUE)  
landcover05 <- as.data.frame(urban2005)
main_urban<-bind_cols(main_urban, landcover05)
urban2006 <- terra::zonal(modis06, test, fun=mean, na.rm=TRUE)  
landcover06 <- as.data.frame(urban2006)
main_urban<-bind_cols(main_urban, landcover06)
urban2007 <- terra::zonal(modis07, test, fun=mean, na.rm=TRUE) 
landcover07 <- as.data.frame(urban2007)
main_urban<-bind_cols(main_urban, landcover07)
urban2008 <- terra::zonal(modis08, test, fun=mean, na.rm=TRUE)  
landcover08 <- as.data.frame(urban2008)
main_urban<-bind_cols(main_urban, landcover08)
urban2009 <- terra::zonal(modis09, test, fun=mean, na.rm=TRUE)  
landcover09 <- as.data.frame(urban2009)
main_urban<-bind_cols(main_urban, landcover09)
urban2010 <- terra::zonal(modis10, test, fun=mean, na.rm=TRUE)  
landcover10 <- as.data.frame(urban2010)
main_urban<-bind_cols(main_urban, landcover10)



main_urban<-main_urban %>% rename(
  name_en=`syr_adm3@data$Name`
)

main_urban<-main_urban %>% rename(
  urban2001= '2001_01_01_LC_Type1')

main_urban<-main_urban %>% rename(
  urban2002= '2002_01_01_LC_Type1'
  )

main_urban<-main_urban %>% rename(
  urban2003= '2003_01_01_LC_Type1'
  )

main_urban<-main_urban %>% rename(
  urban2004= '2004_01_01_LC_Type1'
  )

main_urban<-main_urban %>% rename(
  urban2005= '2005_01_01_LC_Type1')

main_urban<-main_urban %>% rename(
  urban2006= '2006_01_01_LC_Type1')

main_urban<-main_urban %>% rename(
  urban2007= '2007_01_01_LC_Type1')

main_urban<-main_urban %>% rename(
  urban2008= '2008_01_01_LC_Type1')

main_urban<-main_urban %>% rename(
  urban2009= '2009_01_01_LC_Type1')

main_urban<-main_urban %>% rename(
  urban2010= '2010_01_01_LC_Type1')

write.csv(main_urban, "~/ash_trinh_erl_comment_replication/urbanization.csv")

###NOTE: With Completed Data, Appendix 3 Analysis Can Start here#####

abandon_main<-read.dta("~/ash_trinh_erl_comment_replication/main_abandon.dta")


###Table A7: Urbanization, Nighttime Lights and Abandonment
drought.urban.felm <- felm(lights ~ abandon_pct+abandon_pct:urban | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2008  & abandon_main$year <= 2010,], exactDOF=TRUE)
predrought.urban.felm <- felm(lights ~ abandon_pct+abandon_pct:urban | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year < 2008,], exactDOF=TRUE)
all.urban.felm <- felm(lights ~ abandon_pct+abandon_pct:urban | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main, exactDOF=TRUE)

summary(drought.urban.felm)
summary(predrought.urban.felm)
summary(all.urban.felm)


stargazer(all.urban.felm, predrought.urban.felm, drought.urban.felm,
          type='latex',
          title = 'Land Abandonment and Nighttime Lights',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Nighttime Lights',
          covariate.labels = c( 'Land Abandonment Pct.','Land Abandonment Pct. x Urbanization'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)

##### Table A8: Effect of Lagged Nighttime Lights on Abandonment conditional on urbanization####

abandon_main[, lights_lag := shift(lights), by=.(name_en)]

drought.lightsurban.felm.rev <- felm(abandon_pct ~ lights+lights_lag + lights:urban+lights_lag:urban | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2008  & abandon_main$year <= 2010,], exactDOF=TRUE)
predrought.lightsurban.felm.rev <- felm(abandon_pct ~ lights+lights_lag + lights:urban+lights_lag:urban| name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year < 2008,], exactDOF=TRUE)
all.lightsurban.felm.rev  <- felm(abandon_pct~ lights+lights_lag + lights:urban+lights_lag:urban| name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main, exactDOF=TRUE)

summary(drought.lightsurban.felm.rev)
summary(predrought.lightsurban.felm.rev)
summary(all.lightsurban.felm.rev)

nobs(drought.lightsurban.felm.rev)
nobs(predrought.lightsurban.felm.rev)
nobs(all.lightsurban.felm.rev)

stargazer(all.lightsurban.felm.rev, predrought.lightsurban.felm.rev, drought.lightsurban.felm.rev,
          type='latex',
          title = 'Lagged Nighttime Lights and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Land Abandonment Pct.',
          covariate.labels = c( 'Nighttime Lights', 'Nighttime Lights in prior year', 'Nighttime Lights x Urbanization', 'Nighttime Lights in prior year x Urbanization'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)

abandon_main[, abandon_pct_secondyear := shift(abandon_pct, type="lead"), by=.(name_en)]
abandon_main[, abandon_pct_firstyear := shift(abandon_pct_secondyear, type="lead"), by=.(name_en)]

###Table A9: Effect of First Year Land Abandonment on Nighttime Lights conditional on urbanization##

drought.urban.felm_firstyear <- felm(lights ~ abandon_pct_firstyear+abandon_pct_firstyear:urban  | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2006  & abandon_main$year <= 2008,], exactDOF=TRUE)
predrought.urban.felm_firstyear  <- felm(lights ~ abandon_pct_firstyear+abandon_pct_firstyear:urban  | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year < 2006,], exactDOF=TRUE)
all.urban.felm_firstyear  <- felm(lights ~ abandon_pct_firstyear+abandon_pct_firstyear:urban  | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main, exactDOF=TRUE)

summary(drought.urban.felm_firstyear )
summary(predrought.urban.felm_firstyear )
summary(all.urban.felm_firstyear )

nobs(drought.urban.felm_firstyear )
nobs(predrought.urban.felm_firstyear )
nobs(all.urban.felm_firstyear )

stargazer(all.urban.felm_firstyear, predrought.urban.felm_firstyear, drought.urban.felm_firstyear,
          type='latex',
          title = 'Land Abandonment and Nighttime Lights',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Last Year Nighttime Lights',
          covariate.labels = c( 'First Year Abandonment Pct.', 'First Year Abandonment Pct. x Urbanization'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)

###Table A10: Effect of Nighttime Lights on First Year Land Abandonment conditional on urbanization###


drought.lightsurban.felm.rev_firstyear <- felm(abandon_pct_firstyear ~ lights+lights_lag+ lights:urban+lights_lag:urban  | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year >= 2006  & abandon_main$year <= 2008,], exactDOF=TRUE)
predrought.lightsurban.felm.rev_firstyear <- felm(abandon_pct_firstyear ~ lights+lights_lag+ lights:urban+lights_lag:urban  | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main[abandon_main$year < 2006,], exactDOF=TRUE)
all.lightsurban.felm.rev_firstyear  <- felm(abandon_pct_firstyear~ lights+lights_lag+ lights:urban+lights_lag:urban  | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=abandon_main, exactDOF=TRUE)

summary(drought.lightsurban.felm.rev_firstyear)
summary(predrought.lightsurban.felm.rev_firstyear)
summary(all.lightsurban.felm.rev_firstyear)

nobs(drought.lightsurban.felm.rev_firstyear)
nobs(predrought.lightsurban.felm.rev_firstyear)
nobs(all.lightsurban.felm.rev_firstyear)


stargazer(all.lightsurban.felm.rev_firstyear, predrought.lightsurban.felm.rev_firstyear, drought.lightsurban.felm.rev_firstyear,
          type='latex',
          title = 'Lagged Nighttime Lights and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: First Year Abandonment Pct.',
          covariate.labels = c( 'Last Year Nighttime Lights', 'Second Year Nighttime Lights', 'Nighttime Lights x Urbanization', 'Nighttime Lights in prior year x Urbanization'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)

####Interaction of Climatic Stress and Urbanization###

drought.somefelm_abandon_urban <- felm(abandon_pct ~ prcp.dev+temp.dev+prcp.dev:temp.dev+ prcp.dev:urban + temp.dev:urban + prcp.dev:temp.dev:urban | name_en + year  | 0 | admin1_en, data=abandon_main[abandon_main$year >= 2008 & abandon_main$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
predrought.somefelm_abandon_urban <- felm(abandon_pct ~ prcp.dev+temp.dev+prcp.dev:temp.dev+prcp.dev:urban + temp.dev:urban + prcp.dev:temp.dev:urban | name_en + year | 0 | admin1_en, data=abandon_main[abandon_main$year < 2008,], exactDOF=TRUE, psdef=FALSE)
all.somefelm_abandon_urban <- felm(abandon_pct ~ prcp.dev+temp.dev+prcp.dev:temp.dev+prcp.dev:urban + temp.dev:urban + prcp.dev:temp.dev:urban | name_en + year | 0 | admin1_en, data=abandon_main, exactDOF=TRUE, psdef=FALSE)

summary(drought.somefelm_abandon_urban)
summary(predrought.somefelm_abandon_urban)
summary(all.somefelm_abandon_urban)

nobs(drought.somefelm_abandon_urban)
nobs(predrought.somefelm_abandon_urban)
nobs(all.somefelm_abandon_urban)

stargazer(all.somefelm_abandon_urban, predrought.somefelm_abandon_urban, drought.somefelm_abandon_urban,
          type='latex',
          title = 'Precipitation, Temperature, and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Land Abandonment Pct.',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "No", "No", "No")),
          notes = c('Standard errors are clustered on','first-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)

####Figure 3: Marginal Effects Plots of Urban*Nighttime Lights/Urban*Land Abandonment Interactions####

vcov(predrought.urban.felm)
vcov(predrought.lightsurban.felm.rev)

conf=.95
title="Effect of Prior-Year Nighttime Lights"
xlabel="Urban Land Cover (MODIS)"
ylabel="Change in Third-Year Land Abandonment"

# Get coefficients of variables
beta_1 =   -0.0005433
beta_3 = 0.0496943 

# Set range of the moderator variable
# Minimum

min_val = 0
max_val = 1



# Create list of moderator values at which marginal effect is evaluated
x_2 <- seq(from=min_val, to=max_val, by=1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =        .00002946981 + (x_2^2)* .0002839113 + 2*x_2*-.00001081739 

# Standard errors
se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title, xaxt='n',cex.lab=1.4)

axis(side = 1, 
     at = 0:1,# Draw x-axis
     c("None", "Total"))


# Plot estimated effects
lines(y=delta_1, x=x_2)
lines(y=upper_bound, x=x_2, lty=2)
lines(y=lower_bound, x=x_2, lty=2)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)

conf=.95
title="Effect of Third-Year Land Abandonment"
xlabel="Urban Land Cover (MODIS)"
ylabel="Change in Year-of Nighttime Lights"

# Get coefficients of variables
beta_1 =   0.08311
beta_3 =3.56557 

# Set range of the moderator variable
# Minimum

min_val = 0
max_val = 1



# Create list of moderator values at which marginal effect is evaluated
x_2 <- seq(from=min_val, to=max_val, by=1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =         0.1352593 + (x_2^2)* 4.088212 + 2*x_2*-0.3771340 

# Standard errors
se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title, xaxt='n',cex.lab=1.4)

axis(side = 1, 
     at = 0:1,# Draw x-axis
     c("None", "Total"))


# Plot estimated effects
lines(y=delta_1, x=x_2)
lines(y=upper_bound, x=x_2, lty=2)
lines(y=lower_bound, x=x_2, lty=2)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)


###############################################################################
###Appendix 4: Separate Effects on Urban and non-urban population #############
###############################################################################

#### Descriptive Statistics of Population by Urban-non Urban within Agricultural Nawahi###


#setwd("~/ash_trinh_erl_comment_replication/MODIS")
library(terra)
library(raster)
library(sf)

modis01 <- rast('modis01.tif')
m <- c(0, 12.1, 0,  12.9, 13.1, 1,  13.9, 17.1, 0)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
modis01 <- classify(modis01, rclmat)

values(modis01)[values(modis01) < 1] = NA

modis01.bin <- raster('modis01_urban_bin.tif')

modis01<-terra::mask(modis01,fullSYR)

#setwd("~/ash_trinh_erl_comment_replication/SYR_adm")
fullSYR <- vect("SYR_ADM3_noflares.shp")

plot(fullSYR)

rm(modis01.polygons)

#setwd("~/ash_trinh_erl_comment_replication/NewFallow")

fallow2000<-rast("Fallow_2000.tif")
fallow2000<-as.polygons(fallow2000)
plot(fallow2000)

AgricultureSYR <- fullSYR[fallow2000,]


writeVector(AgricultureSYR, "~/ash_trinh_erl_comment_replication/MODIS/AgricultureSYR", filetype="ESRI Shapefile", overwrite=TRUE)

##Clip and Erase commands carried out in ArcGIS on full raster borders
##Replication data includes files, compilation not necessary


SYRAg.urban01<-vect("~/ash_trinh_erl_comment_replication/MODIS/AgricultureSYR_urban.shp")
SYRAg.rural01<-vect("~/ash_trinh_erl_comment_replication/MODIS/AgricultureSYR_rural.shp")


SYRAg.urban01_noflares<-terra::crop(SYRAg.urban01,AgricultureSYR)
SYRAg.rural01_noflares<-terra::crop(SYRAg.rural01,AgricultureSYR)

plot(syr_adm3,col="gray")
plot(AgricultureSYR, col="green", add=TRUE)
plot(SYRAg.urban01_noflares, col="purple", add=TRUE)

#setwd("~/ash_trinh_erl_comment_replication/Nighttime Lights")
lights01<-rast("F142001.v4b_web.stable_lights.avg_vis.tif")
lights02<-rast("F142002.v4b_web.stable_lights.avg_vis.tif")
lights03<-rast("F142003.v4b_web.stable_lights.avg_vis.tif")
lights04<-rast("F162004.v4b_web.stable_lights.avg_vis.tif")
lights05<-rast("F152005.v4b_web.stable_lights.avg_vis.tif")
lights06<-rast("F152006.v4b_web.stable_lights.avg_vis.tif")
lights07<-rast("F152007.v4b_web.stable_lights.avg_vis.tif")
lights08<-rast("F162008.v4b_web.stable_lights.avg_vis.tif")
lights09<-rast("F162009.v4b_web.stable_lights.avg_vis.tif")
lights10<-rast("F182010.v4d_web.stable_lights.avg_vis.tif")

citylights01 <- extract(lights01, SYRAg.urban01_noflares, fun=mean, na.rm=FALSE) 
citylights01 <- as.data.frame(citylights01)
rurallights01 <- extract(lights01, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights01 <- as.data.frame(rurallights01)
citylights02 <- extract(lights02, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights02 <- as.data.frame(citylights02)
rurallights02 <- extract(lights02, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights02 <- as.data.frame(rurallights02)
citylights03 <- extract(lights03, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights03 <- as.data.frame(citylights03)
rurallights03 <- extract(lights03, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights03 <- as.data.frame(rurallights03)
citylights04 <- extract(lights04, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights04 <- as.data.frame(citylights04)
rurallights04 <- extract(lights04, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights04 <- as.data.frame(rurallights04)
citylights05 <- extract(lights05, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights05 <- as.data.frame(citylights05)
rurallights05 <- extract(lights05, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights05 <- as.data.frame(rurallights05)
citylights06 <- extract(lights06, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights06 <- as.data.frame(citylights06)
rurallights06 <- extract(lights06, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights06 <- as.data.frame(rurallights06)
citylights07 <- extract(lights07, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights07 <- as.data.frame(citylights07)
rurallights07 <- extract(lights07, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights07 <- as.data.frame(rurallights07)
citylights08 <- extract(lights08, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights08 <- as.data.frame(citylights08)
rurallights08 <- extract(lights08, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights08 <- as.data.frame(rurallights08)
citylights09 <- extract(lights09, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights09 <- as.data.frame(citylights09)
rurallights09 <- extract(lights09, SYRAg.rural01_noflares, fun=mean,  na.rm=FALSE) 
rurallights09 <- as.data.frame(rurallights09)
citylights10 <- extract(lights10, SYRAg.urban01_noflares, fun=mean,  na.rm=FALSE) 
citylights10 <- as.data.frame(citylights10)
rurallights10 <- extract(lights10, SYRAg.rural01_noflares, fun=mean, na.rm=FALSE) 
rurallights10 <- as.data.frame(rurallights10)


library(dplyr)

urban.names<-as.data.frame(SYRAg.urban01_noflares$Name)
rural.names<-as.data.frame(SYRAg.rural01_noflares$Name)

citylights01<-citylights01 %>% rename(
  citylights2001=F142001.v4b_web.stable_lights.avg_vis
)
rurallights01<-rurallights01 %>% rename(
  rurallights2001=F142001.v4b_web.stable_lights.avg_vis
)
citylights02<-citylights02 %>% rename(
  citylights2002=F142002.v4b_web.stable_lights.avg_vis
)
rurallights02<-rurallights02 %>% rename(
  rurallights2002=F142002.v4b_web.stable_lights.avg_vis
)
citylights03<-citylights03 %>% rename(
  citylights2003=F142003.v4b_web.stable_lights.avg_vis
)
rurallights03<-rurallights03 %>% rename(
  rurallights2003=F142003.v4b_web.stable_lights.avg_vis
)
citylights04<-citylights04 %>% rename(
  citylights2004=F162004.v4b_web.stable_lights.avg_vis
)
rurallights04<-rurallights04 %>% rename(
  rurallights2004=F162004.v4b_web.stable_lights.avg_vis
)
citylights05<-citylights05 %>% rename(
  citylights2005=F152005.v4b_web.stable_lights.avg_vis
)
rurallights05<-rurallights05 %>% rename(
  rurallights2005=F152005.v4b_web.stable_lights.avg_vis
)
citylights06<-citylights06 %>% rename(
  citylights2006=F152006.v4b_web.stable_lights.avg_vis
)
rurallights06<-rurallights06 %>% rename(
  rurallights2006=F152006.v4b_web.stable_lights.avg_vis
)
citylights07<-citylights07 %>% rename(
  citylights2007=F152007.v4b_web.stable_lights.avg_vis
)
rurallights07<-rurallights07 %>% rename(
  rurallights2007=F152007.v4b_web.stable_lights.avg_vis
)
citylights08<-citylights08 %>% rename(
  citylights2008=F162008.v4b_web.stable_lights.avg_vis
)
rurallights08<-rurallights08 %>% rename(
  rurallights2008=F162008.v4b_web.stable_lights.avg_vis
)
citylights09<-citylights09 %>% rename(
  citylights2009=F162009.v4b_web.stable_lights.avg_vis
)
rurallights09<-rurallights09 %>% rename(
  rurallights2009=F162009.v4b_web.stable_lights.avg_vis
)
citylights10<-citylights10 %>% rename(
  citylights2010=F182010.v4d_web.stable_lights.avg_vis
)
rurallights10<-rurallights10 %>% rename(
  rurallights2010=F182010.v4d_web.stable_lights.avg_vis
)
citylights01<-cbind(citylights01,urban.names,2001)
colnames(citylights01)[4]<-"year"
citylights02<-cbind(citylights02,urban.names,2002)
colnames(citylights02)[4]<-"year"
citylights03<-cbind(citylights03,urban.names,2003)
colnames(citylights03)[4]<-"year"
citylights04<-cbind(citylights04,urban.names,2004)
colnames(citylights04)[4]<-"year"
citylights05<-cbind(citylights05,urban.names,2005)
colnames(citylights05)[4]<-"year"
citylights06<-cbind(citylights06,urban.names,2006)
colnames(citylights06)[4]<-"year"
citylights07<-cbind(citylights07,urban.names,2007)
colnames(citylights07)[4]<-"year"
citylights08<-cbind(citylights08,urban.names,2008)
colnames(citylights08)[4]<-"year"
citylights09<-cbind(citylights09,urban.names,2009)
colnames(citylights09)[4]<-"year"
citylights10<-cbind(citylights10,urban.names,2010)
colnames(citylights10)[4]<-"year"

rurallights01<-cbind(rurallights01,rural.names,2001)
colnames(rurallights01)[4]<-"year"
rurallights02<-cbind(rurallights02,rural.names,2002)
colnames(rurallights02)[4]<-"year"
rurallights03<-cbind(rurallights03,rural.names,2003)
colnames(rurallights03)[4]<-"year"
rurallights04<-cbind(rurallights04,rural.names,2004)
colnames(rurallights04)[4]<-"year"
rurallights05<-cbind(rurallights05,rural.names,2005)
colnames(rurallights05)[4]<-"year"
rurallights06<-cbind(rurallights06,rural.names,2006)
colnames(rurallights06)[4]<-"year"
rurallights07<-cbind(rurallights07,rural.names,2007)
colnames(rurallights07)[4]<-"year"
rurallights08<-cbind(rurallights08,rural.names,2008)
colnames(rurallights08)[4]<-"year"
rurallights09<-cbind(rurallights09,rural.names,2009)
colnames(rurallights09)[4]<-"year"
rurallights10<-cbind(rurallights10,rural.names,2010)
colnames(rurallights10)[4]<-"year"

colnames(rurallights01)[3]<-"NAME_EN"
colnames(rurallights02)[3]<-"NAME_EN"
colnames(rurallights03)[3]<-"NAME_EN"
colnames(rurallights04)[3]<-"NAME_EN"
colnames(rurallights05)[3]<-"NAME_EN"
colnames(rurallights06)[3]<-"NAME_EN"
colnames(rurallights07)[3]<-"NAME_EN"
colnames(rurallights08)[3]<-"NAME_EN"
colnames(rurallights09)[3]<-"NAME_EN"
colnames(rurallights10)[3]<-"NAME_EN"

colnames(citylights01)[3]<-"NAME_EN"
colnames(citylights02)[3]<-"NAME_EN"
colnames(citylights03)[3]<-"NAME_EN"
colnames(citylights04)[3]<-"NAME_EN"
colnames(citylights05)[3]<-"NAME_EN"
colnames(citylights06)[3]<-"NAME_EN"
colnames(citylights07)[3]<-"NAME_EN"
colnames(citylights08)[3]<-"NAME_EN"
colnames(citylights09)[3]<-"NAME_EN"
colnames(citylights10)[3]<-"NAME_EN"


colnames(rurallights01)[2]<-"rurallights"
colnames(rurallights02)[2]<-"rurallights"
colnames(rurallights03)[2]<-"rurallights"
colnames(rurallights04)[2]<-"rurallights"
colnames(rurallights05)[2]<-"rurallights"
colnames(rurallights06)[2]<-"rurallights"
colnames(rurallights07)[2]<-"rurallights"
colnames(rurallights08)[2]<-"rurallights"
colnames(rurallights09)[2]<-"rurallights"
colnames(rurallights10)[2]<-"rurallights"

colnames(citylights01)[2]<-"citylights"
colnames(citylights02)[2]<-"citylights"
colnames(citylights03)[2]<-"citylights"
colnames(citylights04)[2]<-"citylights"
colnames(citylights05)[2]<-"citylights"
colnames(citylights06)[2]<-"citylights"
colnames(citylights07)[2]<-"citylights"
colnames(citylights08)[2]<-"citylights"
colnames(citylights09)[2]<-"citylights"
colnames(citylights10)[2]<-"citylights"

citylights<-rbind(citylights01,citylights02,citylights03,citylights04,citylights05,citylights06,citylights07,citylights08,citylights09,citylights10)
rurallights<-rbind(rurallights01,rurallights02,rurallights03,rurallights04,rurallights05,rurallights06,rurallights07,rurallights08,rurallights09,rurallights10)

new.lights<-merge(rurallights,citylights,by=c("NAME_EN","year"))
drops <- c("ID.x","ID.y")
new.lights<-new.lights[ , !(names(new.lights) %in% drops)]
colnames(new.lights)[1]<-"name_en"
setwd("~/ash_trinh_erl_comment_replication")

final_abandon<-merge(abandon_main,new.lights, by=c("name_en","year"))
write.csv(final_abandon,"final_abandon.csv")

###Replication includes files with completed mergin. Compilation not necessary for analysis

###Table A11: Models for Effect of Abandonment on Urban and Rural Nighttime Lights
d.abandon.urban <- felm(citylights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=final_abandon[final_abandon$year >= 2008  & final_abandon$year <= 2010,], exactDOF=TRUE)
d.abandon.rural <- felm(rurallights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=final_abandon[final_abandon$year >= 2008  & final_abandon$year <= 2010,], exactDOF=TRUE)
pd.abandon.urban <- felm(citylights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=final_abandon[final_abandon$year < 2008,], exactDOF=TRUE)
pd.abandon.rural <- felm(rurallights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=final_abandon[final_abandon$year < 2008,], exactDOF=TRUE)
all.urban <- felm(citylights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=final_abandon, exactDOF=TRUE)
all.rural <- felm(rurallights ~ abandon_pct | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin2_en, data=final_abandon, exactDOF=TRUE)

summary(d.abandon.urban)
summary(d.abandon.rural)
summary(pd.abandon.urban)
summary(pd.abandon.rural)
summary(all.urban)
summary(all.rural)



stargazer(all.urban,pd.abandon.urban,d.abandon.urban,all.rural,pd.abandon.rural,d.abandon.rural,
          type='latex',
          title = 'Lights and Land Abandonment',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years - Urban', 'Pre-Drought - Urban', 'Drought - Urban','All Years - Rural', 'Pre-Drought - Rural', 'Drought - Rural'),
          dep.var.caption = 'DV: Average Nighttime Lights',
          covariate.labels = c( 'Land Abandonment'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only nawahi Elklund et. al. (2024)', 'record as having agricultural land',
                    'are included in the model: 184 with rural land and 173 urban land.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)


###Effects of Climatic Stress on Urban Population Change###

drought.felm_full.urban <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=final_abandon[final_abandon$year >= 2008 & final_abandon$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
drought.lm_full.urban <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=final_abandon[final_abandon$year >= 2008 & final_abandon$year <= 2010,])
predrought.felm_full.urban <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=final_abandon[final_abandon$year < 2008,], exactDOF=TRUE, psdef=FALSE)
predrought.lm_full.urban <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=final_abandon[final_abandon$year < 2008,])
all.felm_full.urban <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=final_abandon, exactDOF=TRUE, psdef=FALSE)
all.lm_full.urban <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=final_abandon)

summary(drought.felm_full.urban)
summary(drought.lm_full.urban)
summary(predrought.felm_full.urban)
summary(predrought.lm_full.urban)
summary(all.felm_full.urban)
summary(all.lm_full.urban)


###Effects of Climatic Stress on Rural Population Change###

drought.felm_full.rural <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=final_abandon[final_abandon$year >= 2008 & final_abandon$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
drought.lm_full.rural <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=final_abandon[final_abandon$year >= 2008 & final_abandon$year <= 2010,])
predrought.felm_full.rural <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en  + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=final_abandon[final_abandon$year < 2008,], exactDOF=TRUE, psdef=FALSE)
predrought.lm_full.rural <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=final_abandon[final_abandon$year < 2008,])
all.felm_full.rural <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name_en + as.factor(admin1_en):as.factor(year) | 0 | admin1_en, data=final_abandon, exactDOF=TRUE, psdef=FALSE)
all.lm_full.rural <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=final_abandon)

summary(drought.felm_full.rural)
summary(drought.lm_full.rural)
summary(predrought.felm_full.rural)
summary(predrought.lm_full.rural)
summary(all.felm_full.rural)
summary(all.lm_full.rural)

